home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 06 - 1990 / 06.11 Nov 90 / NeuralNet Estimator⁄EDIT / Ord Least Squares.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-08-03  |  6.3 KB  |  226 lines  |  [TEXT/KAHL]
  1. #include "Neural Network.h"
  2. #include <math.h>
  3.  
  4. extern FILE * Jac;
  5.  
  6. /*-------------------------
  7.     Compute the next step of iteration by solving linear system.
  8. */
  9. OLSbyQRmethod(X,P,D,Y)
  10.     DTypeMatrix * X;    /* pointer to transpose of explanatory data */
  11.     DTypeVector * P;    /* pointer to Pi vector */
  12.     DTypeVector * D;    /* pointer to Diag vector */
  13.     DTypeVector * Y;    /* pointer to vector of dependent variables*/
  14. {
  15.     int sing;
  16.  
  17.     sing = QRDecomposition(X,P,D);
  18.     if(sing==FALSE)
  19.     {    ComputeQY(X,P,Y);
  20.         SolveRbY(X,D,Y);
  21.     }
  22.     return(sing);
  23. }
  24.  
  25.  
  26. /*-------------------------
  27.     Compute Q*Y to get dependent variable for triangularized system, where 
  28.     Q=U(N)*U(N-1)*...*U(1), is the product of N elementary reflecting matrices.  
  29.     Uses output from QR decomposition, see Stewert, p239.
  30. */
  31. ComputeQY(Q,P,Y)
  32.     DTypeMatrix    * Q;    /* matrix output from QR decomposition */
  33.     DTypeVector * P;    /* vector of pi values from QR decomposition */
  34.     DTypeVector    * Y;    /* vector of dependent values */
  35. {
  36.     register int    i;
  37.     register int    k; 
  38.     register int    M; 
  39.     register int    N;
  40.     register DataType    sum;    /* used to temporarily sum terms vector product */
  41.     register DataType    * u;    /* pointer to values used in elementary reflecting matrix */
  42.     register DataType    * y;    /* pointer to array of dependent values */
  43.     register DataType    * pi;    /* pointer to array of pi values from QR decomposition */
  44.     
  45.     HLock(Q->cells);
  46.     HLock(P->cells);
  47.     HLock(Y->cells);
  48.     
  49.     M = Q->cols;
  50.     N = Q->rows;
  51.     pi = *P->cells;
  52.     y = *Y->cells;
  53.  
  54.     for(i=0; i<N; i++, pi++)
  55.     {    u = *Q->cells +i*M;
  56.         sum = 0.0;
  57.         for(k=i; k<M; k++)
  58.             sum += u[k]*y[k];
  59.         sum = sum/(*pi);
  60.         for(k=i; k<M; k++)
  61.             y[k] -= sum*u[k];
  62.     }
  63.     
  64.     HUnlock(Q->cells);
  65.     HUnlock(P->cells);
  66.     HUnlock(Y->cells);
  67. }
  68.  
  69. /*-------------------------
  70.     Algorithm 3.1.3 of Stewart
  71.     Solve linear system Rb=y for b, where R is an upper triangular nonsingular matrix.
  72.     R is stored as its transpose so R(i,j) is at jth row, ith column. 
  73.     The diagonal terms are in seperate vector D as described in Algorithm 5.3.8 of 
  74.     Stewart.  Answer is returned in vector D.
  75. */
  76. SolveRbY(R,D,Y)
  77.     DTypeMatrix    * R;
  78.     DTypeVector    * D;    /* vector of diagonal terms for upper triangular matrix R */
  79.     DTypeVector    * Y;    /* vector of dependent values in linear model */
  80. {
  81.     DataType     * r;
  82.     register int    i;
  83.     register int    j; 
  84.     register int    M; 
  85.     register int    N;
  86.     register DataType * b;    /* pointer to parameter values */
  87.     register DataType * y;    /* pointer to dependent values in linear model */
  88.     register DataType sum;
  89.     register DataType * temp;
  90.     
  91.     HLock(R->cells);
  92.     HLock(D->cells);
  93.     HLock(Y->cells);
  94.     
  95.     r = *R->cells;
  96.     M = R->cols;
  97.     N = D->rows;            /* this also equals R->rows, is the number of parameters */
  98.     b = *D->cells;            /* use the diagonal vector to store parameter estimates */
  99.     y = *Y->cells;
  100.     for(i=N-1; i>-1; i--)
  101.     {    sum = 0.0;
  102.         for(j=i+1, temp = r+i+j*M; j<N; j++, temp = temp + M)
  103.             sum += (*temp)*b[j];
  104.         b[i] = (y[i]-sum)/b[i];
  105.     }
  106.     
  107.     HUnlock(R->cells);
  108.     HUnlock(D->cells);
  109.     HUnlock(Y->cells);
  110. }
  111.  
  112. /*---------------------
  113.     QR decomposition algorithm, see Algorithm 3.8 in 
  114.     Introduction to matrix Computations by G. Stewart, also Algorithm A3.2.1 in
  115.     Numerical Methods for Unconstrained Optimization and Nonlinear Equations by 
  116.     Dennis and Shnabel.  Used to solve linear system Ax=b, where A is (MxM), x is (Nx1).
  117.     
  118.     For coding efficiency the input matrix A is the transpose of the matrix given in
  119.     the statement of the algorithm available in above references.
  120.     
  121.     Assumes more observations than parameters, ie M>N.
  122. */
  123.  
  124. QRDecomposition(A,P,D)
  125.     DTypeMatrix    * A;        
  126.     DTypeVector    * P, * D;
  127. {
  128.     register int k, j, i;
  129.     register int    N;            /* number of rows in A    */
  130.     register int    M;            /* number of columns in A */
  131.     int sing = FALSE;            /* flag for singular A matrix */
  132.     DataType * kth_col;            /* points to row of A, column of original untransposed matrix */
  133.     DataType * pi;                /* pointer to array for the pi values */
  134.     register DataType * diag;    /* pntr to array of diagonal terms for the triangular R matrix */
  135.     register DataType * alpha;    /* temporary pointer to cells in row pointed to by kth_col */
  136.     register DataType * temp;    /* temporary pointer */
  137.     register DataType sign;        /* sign of diagonal term in A matrix */
  138.     register DataType aida, sigma, tau;
  139.     
  140.     HLock(A->cells);
  141.     HLock(P->cells);
  142.     HLock(D->cells);
  143.     
  144.     kth_col  = *A->cells;    /* start off in row zero */
  145.     N = A->rows;        /* number of parameters in linear system, since using transpose */
  146.     M = A->cols;        /* number of observations, since using transpose */
  147.     pi = *P->cells;
  148.     diag = *D->cells;
  149.  
  150.     for(k=0; k<N; k++, kth_col = kth_col + M, pi++, diag++)
  151.     {
  152.         aida = 0.0;                            /* initialize max abs value to zero */
  153.         temp = kth_col+k;
  154.         for(i=k; i<M; i++, temp++ )
  155.         {    tau = fabs(*temp);                /* calculate aida     */
  156.             if(aida < tau) aida = tau;        /* just borrowing tau */
  157.         }
  158.         
  159.         if(aida == 0.0)
  160.         {    *pi = 0.0;        /* column is already triangular */
  161.             *diag = 0.0;    /* see line 2.2T.2 of Alg A3.2.1 from Dennis and Schnabel */
  162.             sing = TRUE;
  163.         }
  164.         else
  165.         {    
  166.             for(i=k, alpha = kth_col + k; i<M; i++, alpha++)
  167.                 *alpha = (*alpha)/aida;
  168.             sigma = 0.0;
  169.             temp = kth_col+k;
  170.             if(*temp>0) sign = 1;                /*         calculate sign term     */
  171.             else sign = -1;
  172.  
  173.             for(i=k; i<M; i++, temp++ )            /*------                        */
  174.                 sigma += (*temp)*(*temp);        /*         calculate sigma            */
  175.             sigma = sign*sqrt(sigma);            /*------                        */            
  176.  
  177.             *(kth_col + k) += sigma;
  178.             *pi = (*(kth_col+k))*sigma;
  179.             *diag = -aida*sigma;
  180.  
  181.             for(j=k+1; j<N; j++)
  182.             {    
  183.                 temp =kth_col+k;                        /*------            */
  184.                 alpha=kth_col+k + (j-k)*M;                /*                    */
  185.                 tau = 0.0;                                /* calculate tau    */
  186.                 for(i=k; i<M; i++, alpha++, temp++)        /*                     */
  187.                     tau += (*temp)*(*alpha);            /*                    */
  188.                 tau = tau/(*pi);                        /*------            */
  189.                 
  190.                 temp =kth_col+k;
  191.                 alpha=kth_col+k + (j-k)*M;
  192.                 for(i=k; i<M; i++, alpha++, temp++)
  193.                     *alpha -= (*temp)*tau;
  194.             }
  195.         }
  196.     }    /*end of for(k=0; k<N; k++, kth_col = kth_col + M, pi++, diag++)*/
  197.     HUnlock(A->cells);
  198.     HUnlock(P->cells);
  199.     HUnlock(D->cells);
  200.     return(sing);
  201. }
  202.  
  203. WriteYXToFile(jac,vector,matrix)
  204.     FILE * jac;
  205.     DTypeVector * vector;
  206.     DTypeMatrix * matrix;
  207. {
  208.     int    i,j;
  209.     DataType * mcell;
  210.     DataType * vcell;
  211.     
  212.     HLock(matrix->cells);
  213.     HLock(vector->cells);
  214.     
  215.     mcell = *matrix->cells;
  216.     vcell = *vector->cells;
  217.     for(j=0; j<matrix->cols; j++)
  218.     {    fprintf(jac,"%.5lf    ",vcell[j]);
  219.         for(i=0; i<matrix->rows; i++)
  220.             fprintf(jac,"%.8lf    ",*(mcell + (i*matrix->cols + j)) );
  221.         fprintf(jac,"\n");
  222.     }
  223.     HUnlock(matrix->cells);
  224.     HUnlock(vector->cells);
  225. }    
  226.